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FOREWORD 


In  this  report  there  are  presented  efficient  methods  for  the  computation  of  velocity 
in  the  wave  train  of  a point  source.  Part  of  the  material  in  the  report  has  been  taken 
from  previous  reports  an-*  part  has  been  developed  in  recent  work.  The  report  is 
intended  to  provide  a systematic  coherent  documentation  for  subroutines.  The 
manuscript  was  completed  by  7 September  197'/. 
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ABSTRACT 

The  velocity  potential  and  the  components  of  velocity  in  the  wave  train  of  a point 
source  are  derived  from  Havelock’s  integral.  Analysis  and  documentation  are  given  for 
the  computation  of  velocity  by  trapezoidal  integration,  by  asymptotic  approximation, 
and  by  integration  by  parts.  Accuracy  and  efficiency  are  least  dependent  upon  the 
depth  of  the  source  in  the  integration  by  parts.  The  three  methods  of  computation 
are  available  for  applications  in  a set  of  subroutines. 


INTRODUCTION 


In  the  computation  of  the  flow  around  a ship,  it  is  assumed  usually  that  the  fluid 
is  incompressible  and  irrotationai.  It  is  assumed  that  the  fluid  is  infinite  in  depth, 
and  that  the  free  sunace  can  be  represented  by  a linear  approximation.  In  one  method 
for  the  computation  of  the  flow  around  a ship,  tne  product  of  source  density  at  the 
surface  of  the  ship  and  the  velocity  of  flow  from  a point  source  is  integrated  over  ttye 
surface  of  the  ship. 

The  classical  approach  to  a surface  disturbance  by  Kelvin1  was  an  application  of 
the  principle  of  stationary  phase  to  a determination  of  the  first  term  of  an  asymptotic 
exoansion  Consideration  of  the  possibility  of  development  into  full  expansions  has 
been  made  for  the  line  of  wave  crests  by  Peters*  and  by  Ursell3.  Development  of  the 
full  asymptotic  approximation  has  been  completed  in  the  meantime >T.  The  utilization 
of  the  asymptotic  approximation  is  limited  to  the  far  field. 

Fourier  integrals  for  the  velocity  potential  in  the  wave  train  of  a point  source  have 
been  derived  for  infinite  depth  by  Havelock*  and  for  finite  depth  by  Lunde*.  The 
components  of  velocity  are  known  also  as  the  derivatives  of  Green's  function  as  defined 
by  Wehausen*. 

The  velocity  in  the  flow  from  a point  source  is  expressed  as  the  sum  of  three  terms. 
The  first  term  is  the  flow  from  the  source  in  an  unbounded  fluid,  the  second  term  is 
the  flow  from  an  image  source  above  the  free  surface,  and  the  third  term  is  the  flow 
in  the  wave  train.  The  third  term  often  is  expressed  as  the  sum  of  a single  integral 
which  represents  a free  wave,  and  a double  integral  which  represents  a bound  wave. 
The  single  integral  by  itself  satisfies  the  linearized  free-boundary  conditions,  while 
the  double  integral  combines  with  the  first  two  terms  to  satisfy  the  linearized 
free-boundary  conditions.  The  amount  by  which  the  single  integral  should  be  added 
can  be  determined  by  an  analysis  of  the  transient  wave  pattern  of  a source  which 
starts  at  zero  time  and  moves  thereafter  at  a constant  rate. 

Components  of  velocity  in  physical  space  are  derived  from  Fourier  integration  in 
wave  number  space.  Either  Cartesian  coordinates  or  polar  coordinates  may  be  used 
in  either  space  for  the  expression  of  functions.  The  usual  choice  is  Cartesian  coordinates 
in  physical  space  and  polar  coordinates  in  wave  number  space. 

The  single  integral  is  just  the  contribution  from  integration  on  a small  circle  at  a 
singularity  in  the  complex  plane  of  the  radial  coordinate  The  double  integral  is  the 
Cauchy  principal  value  of  an  integral  whose  integrand  is  infinite  at  the  singularity. 
The  single  integral  and  the  double  integral  can  be  combined  into  one  complex  integral 
in  the  complex  plane  of  the  radial  coordinate.  The  integrand  of  the  complex  integral 
is  analytic  along  the  path  of  integration.  It  is  necessary  that  the  integrand  be  analytic, 
because  the  integral  is  the  steady  state  limit  of  the  transient  wave  which  has  an 
analytic  integrand.  It  is  not  necessary  to  resort  to  ficticious  frictions  which  have  no 
physical  meaning  in  an  inviscid  fluid. 

Wherever  the  integrand  is  analytic,  then  in  accordance  with  the  Cauchy  theorem, 
the  path  of  integration  may  be  deformed  into  any  other  path  with  the  same  limits  of 
integration.  The  path  of  integration  with  respect  to  either  coordinate  may  be  moved 
out  into  the  complex  plane  for  that  coordinate.  Complex  paths  have  been  proposed 
several  times  in  the  literature7'*.  Many  years  ago  it  was  recognized  that  radial 
integration  in  wave  number  space  is  equivalent  to  the  evaluation  of  the  complex 
exponential  integral  for  which  we  have  standard  methods  of  computation. 
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Even  after  the  radial  integration  in  wave  number  space,  there  remains  an  azimuthal 
integration  in  wave  number  space.  The  integrand  is  cyclical  with  respect  to  azimuth 
angle,  and  the  trapezoidal  rule  is  the  high  accuracy  rule  for  numerical  quadrature. 
The  integrand  is  partly  monotonic  and  partly  oscillatory,  and  a large  number  of 
intervals  is  required  by  the  trapezoidal  rule,  unless  the  integration  is  for  a position 
in  the  neighborhood  of  the  source. 

A breakthrough  was  achieved  when  the  quadrature  was  replaced  with  an  integration 
by  parts  right  through  any  number  of  fluctuations  of  the  integrand.  Convergent, 
approximation  is  used  in  connection  with  the  integration  by  parts.  The  convergent 
approximation  is  useful  for  both  the  near  field  and  the  far  field. 

The  evaluation  of  the  Green  function  and  its  derivatives  has  received  enough 
attention  in  the  literature  to  justify  the  release  of  formal  subroutines  for  the 
computation  of  velocity  potential  and  the  components  of  velocity  in  the  wave  train 
of  a point  source.  Such  Subroutines  have  been  under  development  in  this  laboratory14'17 
for  a number  of  years.  The  present  report  is  a concise  documentation  of  the  assumptions 
aid  the  formulations  which  are  basic  to  the  subroutines  Important  parts  of  the 
documentation  are  data  on  accuracy  and  efficiency 

The  specification  of  accuracy  for  a subroutine  depends  upon  whether  the  subroutine 
is  intended  for  a specific  application  or  for  general  utilization.  The  level  of  accuracy 
in  the  present  project  is  eight  decimal  digits,  while  the  computer  essentially  is  a CDC 
6600.  In  single  precision  this  computer  operates  on  48- bit  numbers.  Most  computers 
operate  on  numbers  with  a fixed  number,  of  binary  bits.  A gauge  of  accuracy  for 
arithmetic  operations  is  one  unit  in  the  lowest  order  bit.  Larger  errors  accumulate 
during  the  computation  of  a function.  There  is  inherent  error  which  arises  from  error 
in  the  argument  and  there  are  truncation  ei  rnrs  and  rounding  errors  which  arise 
from  inaccuracy  of  computation.  A sensible  policy  is  to  keep  the  truncation  errors 
and  the  rounding  errors  to  the  same  level  as  the  inherent  errors. 

The  absolute  accuracy  of  a variable  is  expressed  by  the  actual  value  of  the  error 
in  the  variable,  while  the  relative  accuracy  of  the  variable  is  expressed  by  the  ratio 
between  the  valua  of  the  error  and  the  value  of  the  variable.  Whether  the  specification 
of  accuracy  can  be  a uniform  absolute  error  or  a uniform  relative  error  depends  upon 
the  nature  of  the  function.  A uniform  absolute  error  is  appropriate  for  a periodic 
function,  but  a uniform  relative  error  is  appropriate  ior  a monotonic  function.  For 
more  complicated  functions  the  absolute  error  may  be  a constant  fraction  of  a 
reference  function  which  matches  the  computed  unction  only  at  a limited  number 
of  maxima. 

The  number  of  cycles  of  a computation  loop  influences  the  accuracy  of  computation. 
The  loop  may  be  terminated  when  a sum  becomes  constant,  when  a preset  tolerance 
is  met,  or  when  a preset  number  of  cycles  have  occurred.  Termination  for  a constant 
sum  gives  the  full  accuracy  of  the  computer  but  may  require  an  expensive  sensing 
operation.  More  economy  is  possible  if  a simple  empirical  formula  can  be  used  to 
compute  the  number  of  cycles,  but  the  derivation  of  an  empirical  formula  may  require 
a program  of  test  runs  on  a computer. 

In  the  computation  of  the  components  of  velocity  the  gauge  of  accuracy  is  the 
square  root  of  the  sum  of  squares  of  the  errors  in  the  three  components  of  velocity. 
At  small  orders  or  arguments  the  error  diminishes  smoothly  with  order  or  argument 
and  is  dominated  by  truncation  error.  At  large  orders  or  arguments  the  error  fluctuates 
with  an  amplitude  Independent  of  order  or  argument  and  is  dominated  by  rounding 
error.  Transitions  between  modes  of  computation  are  designed  so  as  to  balance 
truncation  error  against  rounding  error. 
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Empirical  curves  have  been  constructed  to  express  the  error  and  the  time  for  each 
of  three  methods  of  computation  The  curves  can  indicate  only  trends  in  error  or 
time  because  they  are  based  upon  actual  runs  on  the  CDC  6600  computer  with  time 
sharing  and  optimization  The  error  is  subject  to  statistical  fluctuations  and  the  time 
depends  upon  the  presence  of  other  jobs  in  the  computer.  Under  the  conditions  of 
the  actual  runs  the  cost  of  computation  was  17^  per  second,  but  in  other  installations 
the  cost  would  differ  in  accordance  with  costing  algorithms. 


CARTESIAN  COORDINATES 

Various  authors  on  ship  hydrodynamics  have  used  various  coordinate  systems.  The 
choice  of  coordinates  for  the  present  analysis  is  the  standard  coordinate  system  in 
aeronautics  and  oceanography.  The  forward  direction  is  the  direction  of  motion  of  the 
source  through  still  fluid  The  origin  of  coordinates  is  above  the  source  at  the 
undisturbed  surface  of  the  fluid.  The  x-coordinate  is  the  distance  forward,  the 
y-coordinate  is  the  distance  to  the  right,  and  the  x-coordinate  is  the  distance 
downward.  The  arrangement  of  the  Cartesian  coordinates  is  illustrated  in  Figure  1 


FREE  SURFACE 

In  accordance  with  the  usual  assumption  of  incompressible,  irrotational  flow,  the 
velocity  in  the  fluid  is  the  negative  gradient  of  a velocity  potential  which  satisfies 
Laplace's  equation.  Under  steady  state  conditions  the  potential  in  the  moving  reference 
frame  is  given  by  the  expression 


Uz  + <p{x.  y.  x)  (1) 

where  Vz  is  the  potential  for  uniform  flow  in  a direction  opposite  to  the  motion  of 
the  source  and  ip(x.  y,  x)  is  the  potential  for  the  local  disturbance  in  the  wave  train. 
The  Cartesian  components  u.  v,  iv  of  local  velocity  are  given  by  the  equations 


•u  =* 


.dx 


Let  the  configuration  of  the  free  surface  be  expressed  by  the  equation 

* + <T(*.  y)  = o 


(2) 

(3) 


For  steady-state  conditions  the  velocity  at  the  free  surface  is  tangent  to  the  surface 
and  the  potential  obeys  the  boundary  equation 


dip  d{  dip 
dy  dy  dz 


Neglect  of  terms  of  second  order  leads  tc  the  boundary  equation 


+ 0 
dx  dz 


(4) 


(5) 
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The  motion  of  the  fluid  is  determined  by  the  Bernoulli  equation 

+ — - gz  = | U* 

P 

where  p is  the  density,  p is  the  pressure,  and  g is  the  acceleration  of  gravity.  At  the 
free  surface  the  pressure  is  constant.  Neglect  of  terms  of  second  order  leads  to  the 
boundary  equation 

+ 9(  - o (7) 


The  free  surface  is  eliminated  from  Equations  (5)  and  (7)  by  differentiation  to  give 
the  equation 


6*9  dip 

Jx*  ~ “"dz 


0 


(8) 


where  ca  is  a critical  wave  number  and  is  defined  by  the  equation 


*o 


9_ 

U* 


(9) 


Along  the  centerline  behind  the  source  the  wave  length  A of  the  transverse  waves  is 
given  by  the  equation 


(10) 


Inasmuch  as  the  velocity  field  is  symmetric  with  respect  to  the  x-axis,  it  is  sufficient 
to  '.;—***  the  analysis  to  positive  values  of  y 


FOURIER  ANALYSIS 

The  computation  is  founded  on  the  assumption  that  velocity  can  be  expressed  by 
a Fourier  transform  The  two-dimensional  Fourier  transform  is  given  by  the  equations 


A(o.  0)  - ~ J J F(x.  y)  *-«•*♦**>  dx  dy 

F(x.  y)  = JJ  A(a.  0)  e4*"**’  da  d0 


(11) 


(12) 


where  x.y  are  Cartesian  coordinates  in  physical  space,  a,  0 are  Cartesian  coordinates 
in  wave  number  space.  F(x.y)  is  a function  in  physical  space,  and  A(a.0)  is  the  Fourier 
amplitude  in  wave  number  space  A potential  p(z,  y,  z)  of  a source  distribution  on  the 
plane  z = h.  can  be  constructed  with  the  equation 

?(z.  y.  z)  = JJa(o.  f)  ,-Jtt*-***—*,*  da  (13) 

which  gives  a solution  of  Laplace's  equation  wherever  z * h.  The  vertical  component 


1 


J 


i 


4 


of  velocity  at  x - h is  given  by  the  equation 


Oz 


j J n/o*  +■  fi*  A(a,  p)  «'<—*>  da  dp  (s-h)  (14) 


where  the  sign  is  plus  if  * > h and  the  sign  is  minus  if  z < h.  In  accordance  with  the 
Gauss  theorem  the  difference  in  - d<p/Pz  on  opposite  sides  of  the  plane  z = h is  4ff<7(x,  y) 
where  o(z,  y)  is  the  surface  source  density  at  the  plane.  Thus  the  Fourier  amplitude 
for  an  arbitrary  distribution  of  source  density  is  given  by  the  equation 


') ' J J *'*■ y>  * Jv 


(15) 


For  a unit  source  at  the  origin  the  Fourier  amplitude  is  given  just  by  the  equation 

1 


A(a,p) 


2n\/a*  + fi* 


(16) 


The  amplitude  varies  inversely  as  the  radial  distance  in  wave  number  space. 

Let  k.  6 be  polar  coordinates  in  wave  number  spuce.  The  coordinates  are  related  by 
the  equations 


a = < cos  0 p - k sin  0 

The  two-dimensional  Fourier  transform  is  given  by  the  equations 


dx  dy 

F(x,  y)  - J ^A(k.  8)  *<«<—• •*v.'»«>  KdKdP 


For  the  unit  source  at  the  origin  the  Fourier  amplitude  is  given  by  the  equation 

1 


A(k.  0) 


2tt« 


(17) 

(18) 

(19) 

(20) 


The  potential  (p,(x.y.z)  for  the  isolated  unit  source  is  given  by  the  equation 

y.  z)  = ~ j J ^ dg  (21) 

The  potential  of  the  point  source  does  not  satisfy  the  free-boundary  conditions.  An 
additional  potential  is  added  in  order  to  meet  the  free-boundary  conditions.  The 
Fourier  amplitude  of  the  additional  potential  is  given  by  the  equation 


A(k,$) 


1 (e0  + k cos ZB) 


P.7TK  (<„  - K COS *9) 

The  natural  expansion  into  partial  fractions  is  expressed  by  the  equation 

1 


A(k.  6)  » - - — + 


2tt<c  ttk(k0  - k cos *0) 


(22) 


(23) 


The  first  term  of  the  expansion  is  the  amplitude  of  a negative  image  source  The 
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potential  y,  z)  for  the  image  source  is  given  by  the  equation 


*«(*.  V.  *) 


jl  r r 
2tr  J.ff  Jf 


Thus  the  potential  ? of  the  surface  disturbance  is  given  by  the  equation 

*{x.  y.  *)  = *,(*,  y.  z)  + pt(x.  y.  x)  -t-  ?j(x,  y,  *)  (25) 

where  ?t  is  the  potential  of  the  source  in  an  unbounded  fluid.  <pt  is  the  potential  of 
an  image  source  over  the  free  surface,  and  ip?  is  the  potential  of  the  surface  wave 
train. 

The  potential  <pt  is  given  by  the  equation 

f (26) 

|z*  + y*  ♦ (z  - fi)*)1 

Its  derivatives  are  given  by  the  equations 


0<Pi  _ 

X 

dx 

|x*  + y*  + (z  - h)*)* 

\*f) 

1 

V 

(oa\ 

ay  = 

\xz  + y1  + (z  - h)z\* 

dipt 

z - h. 

/OQ\ 

dx  ” 

(x*  4-  y*  + (*  - h)*)* 

which  are  used  explicity  in  the  computation  of  flow  from  the  point  source. 
The  potential  <pt  is  given  by  the  equation 


<Pt  i 

|x*  + v* +(*.*)*)* 


Its  derivatives  are  given  by  the  equations 


jx*  + y*  + <z  + h)*j* 

* 1 02) 

dip,  z + h.  . 

— = 5 (33) 

91  |z*  + :y*  + <z  + A)*|* 

which  are  used  explicitly  in  the  computation  of  flow  from  the  image  source 

The  Fourier  analysis  of  the  steady  state  gives  only  the  Cauchy  principal  value  of  a 
double  integral.  The  full  expression  for  the  potential  of  the  steady  wave  is  derived  in 
Appendix  A f^om  the  limit  which  is  approached  by  a transient  wave. 


The  potential  <ft  is  given  by  the  equation 


. f+,r  cos  0 

P+w  rm  o ~«{*+JO+i*(*co*#+v*in#) 

- I 

**  J-w  Jo  *0  * 


*0 

2~(*+l»)+i (xeo»#4-y*in#) 

oo*  ‘"d  go*2# 


<ifl 


f -PV- 


cos  *0 


dK  d£ 


-■*  Vo  ''■0 

The  integration  can  be  simplified  through  the  substitutions 

*o 


<5  = 


and 


cos*0 


*0 


(z  + h)  - i(x  cos  0 + y sin  0)  f 


u 


cos*©  (z  + h)  - i(x  cos  8 + y sin  0) 


(34) 


(35) 


(36) 


where  u is  a new  Variable  of  integration.  The  function  e-a  itself  satisfies  both  Laplace's 
equation  and  the  free-surfaee  boundary  equation.  It  is  added  in  just  the  right  amount 
to  make  the  wave  train  trail  behind  the  source  if  the  path  of  integration  with  respect 
to  u proceeds  along  a straight  line  in  the  complex  plane  from  toward  the 

origin,  bypasses  the  origin  on  a half  circle  of  smnll  radius  at  the  origin,  and  continues 
on  a straight  line  to  u * <5.  Whether  the  half  circle  is  on  the  right  or  on  the  left  of 
the  origin  depends  upon  the  sign  of  cos  0.  The  path  of  integration  in  each  case  is 
illustrated  in  Figure  2.  The  potential  of  the  surface  wave  is  given  by  the  equation 

?.<*•  ».*>-*?  J__  jjpjj  ~ *•  ■»  <37> 

Differentiation  of  this  potential  requires  the  equation 


which  may  be  integrated  by  parts  to  give  the  equation 


e . / « 

du  ■-«*•(  — = du 

u J—'u 


The  derivatives  of  ifij  are  given  by  the  equations 


(39) 


(40) 

(41) 

(42) 


When  0 is  replaced  in  the  integrands  by  0 ± n the  integrands  are  replaced  by  their 
complex  conjugates.  The  real  parts  of  integrals  from  -tt  to  +u  are  equal  to  twice  the 
real  parts  of  integrals  from  - to  + Jit. 
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It  is  convenient  to  introduce  new  parameter*  n,  u which  are  deftned  by  the  equations 
H m - (*  ♦ h)  u » * co*  0 * y »in  0 (43) 

whence  the  parameter  6 is  given  by  the  equation 

6 ’ ~ {ti  * ^ (44) 

The  parameter  8 has  always  a positive  real  part,  but  it  may  have  a positive  or  negative 
imaginary  part  according  to  the  value  of  8 

It  is  convenient  to  introduce  a new  parameter  t which  is  defined  by  the  equation 

t - tan  8 (45) 

The  potential  is  given  by  the  equetion 


and  the  derivatives  of  are  given  by  the  equation* 


J"  ^777,  e-sj 

'*  e“ 

— j du  dt 
— u* 

(47) 

,a 

— I d-u.  dt 

— u 

(48) 

j * */.>•■> --M 

gU 

— j du  dt 
— « 

(49) 

■ while  the  parameter  & is  given  by  the  equation 

1 a * - «om(  i + <*)  - «o(*  + yf )■ 

v/l  + f* 

(50) 

If  the  coordinates  h.  x,  y,  z are  scaled  through  multiplication  by  /c„  before  computation, 
then  the  analysis  may  be  simplified  to  the  case  for  which  = 1. 


COMPLEX  PHASE 

Stationary  Phase 

The  complex  phase  8 is  given  in  terms  of  f by  the  equation 

6 = - n(\  * f*)  - t(x  + yf)>/i  + <*  (51) 

There  are  two  branches  of  the  Riemann  surface  over  the  t- plane  The  primary  branch 
contains  the  real  axis,  on  which  the  radical  is  positive  The  secondary  branch  is  reached 
along  any  path  which  encircles  one  of  the  branch  points  at  ±t.  All  members  of  the 
equation  may  be  brought  to  the  same  side,  and  the  expression  so  obtained  may  be 
cleared  of  the  radical  through  multiplication  by  the  same  expression  with  the  sign  of 
the  radical  reversed  The  result  is  the  quartic  equation 

x*  + (m  + 6)*  + Zxyt  + (x*  + y*  + 2m(m  + <5)|f*  + 2xyfJ  + (y*  -*•  nl)t*  = 0 (52) 

There  are  four  branches  of  the  Riemann  surface  over  the  6- plane, 
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j 


f 


l 


The  complex  phase  i is  stationary  where  f satisfies 
— o , (V  + rt  4 2yt ») 

«U  " M vT^Ta 


the  equation 
= 0 


(53) 


The  terms  of  the  equation  may  be  transposed  and  squared  to  clear  the  equation  of 
the  radical.  The  result  is  the  quartic  equation 

V*  + 2 xyt  + (**  + 4y*  + 4//.*)f*  + 4xyf*  + 4(yJ  + m*)<4  - 0 (54) 

which  has  real  coefficients  and  has  two  complex  conjugate  pairs  of  roots.  One  pair  of 
roots  is  associated  with  the  primary  branch  and  the  other  puir  of  roots  is  associated 
with  the  secondary  branch.  The  roots  of  the  quartic  equation  could  be  found  by  the 
Ferrari  method,  but  better  control  over  identification  and  accuracy  is  provided  by 
Newton -Raphson  iteration  with  the  irrational  equation 

The  connection  between  the  complex  d-plane  and  the  complex  f-plane  is  illustrated 
by  Figures  3 and  4,  where  solid  lines  are  various  branches  of  the  real  axis  in  the 
f-plane,  and  dots  and  asterisks  mark  the  points  of  stationary  complex  phase.  In  the 
d-plane.  the  upper  curve  belongs  to  the  primary  branch,  while  the  lower  curve  belongs 
to  the  secondary  branch.  In  the  f-plane.  the  two  lines  which  cross  the  imaginary  axis 
belong  to  the  primary  branch,  while  the  two  lines  which  bend  back  from  the  imaginary 
axis  belong  to  the  secondary  branch. 


c t 


Initial  Approximation 

Solution  of  the  Equation  dd/df  = 0 is  represented  in  an  Argand  diagram  by  the 
vanishing  of  the  sum  of  the  three  vectors 

- ZuitVlTt*  xt  Zy(\  4-  f*)  (55) 

These  vectors  cannot  have  a zero  sum  for  a real  value  of  f unless  x and  y both  are 
zero.  The  first  and  second  vectors  cannot  be  collinear  unless  f is  imaginary  with 
absolute  value  greater  than  1,  in  which  care  the  third  vector  cannot  be  collinear  with 
the  other  two.  The  three  vectors  can  have  a vanishing  combination  only  if  the  first 
and  third  have  components  of  opposite  sign  in  the  direction  perpendicular  to  the 
second  vector  The  limit  of  possibility  is  where  the  second  and  third  vectors  are 
collinear,  in  which  case  f*  and  | form  two  sides  of  a triangle  while  the  third  side  lies 
in  the  direction  of  f.  The  angle  which  f*  makes  with  the  real  axis  is  equal  to  twice 
that  angle  which  t makes,  either  with  the  real  axis,  or  with  the  side  The  triangle 
has  two  angles  equal  and  isjsosceles.  The  equal  sides  are  f*  and  g.  whence  the  limiting 
absolute  value  of  f is  \/'Jz.  . 

If  the  absolute  value  of  t is  less  than  1/vlj,  then  V 1 + tz  deviates  from  the  direction 
of  the  real  axis  by  at  most  gw  The  first  vector  (with  negative  fi)  lies  to  the  left  of 
the  direction  of  t,  whereas  the  third  vector  (with  positive  y)  lies  on  the  positive  real 
side  of  the  direction  of  f The  first  and  third  vectors  oppose  each  other  only  if  the 
imaginary  part  of  t is  positive 

If  the  absolute  value  of  f is  more  than  1/V2,  then  Vl  + f*  still  lies  on  the  positive 
real  side  of  the  imaginary  axis,  provided  the  f-plane  :s  cut  outward  along  the  imaginaiy 
axis  from  ±i.  The  first  vector  still  lies  to  the  left  of  the  direction  of  f,  but  the  third 
vector  now  lies  on  the  negative  real  side  of  the  direction  of  f.  The  first  and  third 
vectors  oppose  each  other  only  if  the  imaginary  part  of  f is  negative. 
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It  may  be  concluded  that  the  equation  di/dt  = 0 has  only  two  pnmary  roots,  of 
which  one  lies  on  the  positive  imaginary  side  of  the  real  axis  inside  a circle  of  radiv« 
1/Vz  while  the  other  lies  on  the  negative  imaginary  side  of  the  real  axis  outside  a 
circle  of  radius  i/v2 

In  the  limit  when  ^<0  the  equation  d6/dt  * 0 is  simplified  to  the  equation 

{ + & + <*-°  (56) 

The  rsots  of  this  equation  are  given  by  the  equations 

- x - v'x*  - 8y*  , - * * 'fx*  - By* 

f, — ft (57) 

Tb,  roots  lie  on  the  real  axis  when  |x|  > N/8|y|.  while  the  roots  lie  on  the  circle  of 
r.tdius  l/v'z  when  !x|  < \/8'y|. 

In  the  limit  when  x ■ 0 the  quartic  equation  is  simplified  to  the  equation 


My*  + M*) 


+ t*  + t*  = 0 


The  roots  of  this  equation  are  given  by  the  equations 


'‘"*75  [“^*771]' 


i r.  m 11 

Vzl'~  v^T?J 


or  by  their  complex  conjugates.  They  are  also  the  roots  of  the  quadratic  equation 


♦ t + f*  = o 

i V Vy1 + |i‘ 


which  is  derived  from  the  polynomial  (f  - t,)(t  - ft). 

In  the  limit  when  y*0  the  quartic  equation  is  simplified  to  the  equation 

4m* 


The  roots  of  this  equation  are  given  by  the  equations 


n/x*  + 4m* 


f»  * * - — 


or  by  their  complex  conjugates  They  are  also  the  roots  of  the  quadratic  equation 

* 

- Vx*  * 471*  f + f*  = 0 (6: 

which  is  derived  from  the  polynomial  (t  - l,)(f  - f*) 

In  the  limit  when  t -*  0,  th3  equation  d6/dt  = 0 is  simplified  to  the  equation 


which  is  derived  from  Equation  (53)  after  neglect  of  the  radical. 
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In  the  limit  when  the  equation  d6/dt  » 0 is  simplified  to  the  equation 

which  is  derived  from  Equation  (53)  after  the  radical  has  been  replaced  by  the  limiting 
approximation. 

vTT<«  ~ * <(  1 f (66) 

The  sign  is  * according  to  whether  the  sign  of  z is  ± 

The  three  special  Equations  (56).  (60).  (63)  may  be  derived  from  the  general  equation 

+ 1 j77  v ]<  + t,.0  (67) 

v y*  ♦ n*  L V*  + M ’ y/y*  + (i‘  J 

which  is  also  compatible  with  the  two  limiting  Equations  (64),  (65)  The  quadratic 
Equation  (67)  is  synthetic,  but  it  gives  the  correct  roots  on  the  real  axis,  on  the 
imaginary  axis,  on  the  circle  of  radius  1/v  2,  and  it  gives  good  enough  roots  elsewhere 
to  serve  as  a good  starting  point  for  the  Newton -Raphson  iteration  of  the  irrational 
Equation  (53) 


TRAPEZOIDAL  INTEGRATION 

In  the  straightforward  evaluation  of  the  Fourier  integrals,  the  radial  integration  is 
expressed  by  the  equation 


■ r * 

j—  u 


for  which  Ft(6)  is  obtained  by  reference  to  a subroutine  for  the  complex  exponential 
integral  The  subroutine  gives  the  integral  along  a path  which  does  not  cross  the 
positive  real  axis.  A correction  must  be  applied  when  the  path  of  integration  does 
cross  the  positive  real  axis.  The  correction 

+ 2-ne'4  (69) 

is  applied  when  fm  6 > 0 The  radial  integration  is  performed  at  equal  intervals  of  the 
azimuth  angle,  and  then  the  azimuthal  integration  is  completed  by  application  of  a 
rule  of  quadrature 

In  order  to  obtain  information  about  amplitude  and  phase,  the  quadrature  is  applied 
through  half  a cycle  of  the  azimuth  angle  The  real  part  is  doubled  and  the  imaginary 
part  is  cancelled  when  the  quadrature  is  continued  through  a full  cycle  of  the  azimuth 
angle  Different  quadrature  rules  are  appropriate  for  the  real  part  and  the  imaginary 
part. 

A function  which  is  continuous  and  periodic  can  be  expressed  in  terms  of  its 
argument  by  a Fourier  series  The  coefficients  in  an  infinite  series  are  derived  from 
integration  of  the  products  of  the  function  and  the  sines  or  cosines  of  multiples  of 


the  argument.  The  coefficients  in  a finite  series  are  derived  from  summations  of  the 
products  of  the  function  and  the  sines  or  cosines  of  multiples  of  equally  spaced  values 
of  the  argument  Orthogonality  of  the  trigonometric  terms  is  true  for  all  orders  in 
integration,  but  is  not  true  for  all  orders  in  summation  If  there  are  ,V  discrete  values 
of  the  argument  then  the  summation  of  the  products  of  sines  or  cosines  is  nonzero 
for  orders  whose  sum  is  a multiple  of  N Thus  computed  values  of  the  coefficients  for 
the  finite  senes  differ  from  the  true  values  of  coefficients  for  the  infinite  series  as 
the  result  of  an  aliasing  of  the  trigonometric  terms  The  number  of  coefficients  in  the 
finite  series  must  be  equal  to  the  number  of  discrete  arguments,  and  the  terms  of 
the  finite  series  with  sines  and  cosines  are  of  order  not  higher  than  JA'  The  aliasing 
error  in  the  constant  of  the  finite  series  is  of  order  at  least  as  high  as  A 

All  terms  in  either  senes  except  the  constants  vanish  identically  when  either  series 
is  integrated  through  one  cycle  The  constant  in  the  finite  series  is  derived  from  the 
sum  of  the  values  of  the  function  at  the  equally  spaced  values  of  the  argument  Thus 
the  rule  for  Fourier  integration  with  respect  to  a cyclical  variable  is  just  the  trapezoidal 
rule  An  assessment  of  the  accuracy  of  the  trapezoidal  rule  requires  an  analysis  of 
the  aliasing  errors  In  the  present  case  of  Fourier  integration  the  coefficients  of  the 
real  part  dimmish  almost  exponentially  with  increasing  order  The  error  of  integration 
is  determined  almost  entirely  by  the  aliased  term  of  lowest  order  Because  of  the 
rapid  decrease  in  the  coefficients  with  increasing  order  the  trapezoidal  rule  is  the 
high  accuracy  rule  for  the  real  part.  In  the  present  case  the  coefficients  of  the 
imaginary  part  diminish  with  order  almost  in  accordance  with  a quadratic  polynomial 
in  the  reciprocal  of  the  order  The  aliasing  error  is  half  as  large  and  reversed  in  sign 
when  the  values  of  the  argument  are  shifted  to  midpoints  between  arguments  Inasmuch 
as  the  Simpson  rule  is  equivalent  to  a superposition  of  two  trapezoidal  rules  with 
shifted  arguments,  the  aliasing  errors  almost  cancel.  The  Simpson  rule  is  a high 
accuracy  rule  for  the  imaginary  part. 

The  accuracy  of  the  integrations  for  the  real  part  was  confirmed  by  comparative 
computations  in  which  the  number  of  intervals  of  trapezoidal  integration  was  increased 
in  steps  until  the  error  was  reduced  to  rounding  error.  The  determination  of  error 
was  repeated  over  enough  coordinates  in  physical  space  to  indicate  the  dependence 
of  error  upon  position  An  empirical  formula  for  the  estimation  of  the  number  of 
intervals  for  eight-digit  accuracy  is  given  by  the  equation 


27  Vx*  - u* 

f + r-S M7  - 54<?  4 ?<7l)  (70) 

Ml 

where  the  argument  q is  given  by  the  equation 

7= 1-  (71) 

**  + !/* 

Inasmuch  as  the  formula  has  only  four  variable  terms,  it  can  be  made  to  satisfy  a 
1 specification  of  accuracy  at  only  four  points  Almost  everywhere  else  the  level  of 

\ ’ accuracy  should  be  better  than  the  specified  level 

I The  dependence  of  the  error  and  the  time  on  position  is  illustrated  by  dashed 

! | curves  in  Figures  5 and  6 Although  the  trapezoidal  integration  is  not  useful  at  small 

I . depth,  it  does  provide  data  for  checkouts  at  moderate  depth 
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ASYMPTOTIC  APPROXIMATION 

Exponential  Integral 

If  |<5|  is  targe  throughout  the  integration  with  respect  to  ®.  then  the  exponential 
integral  ia  given  by  the  asymptotic  approximation 

«'4  f — (/nt iSO)  (72) 
J—  u »-•  0 

when  the  path  of  integration  does  not  encircle  the  origin,  but  it  is  given  by  the 
approximation 

t~*  f*  — du~  £ “j^+2iris"*  (Jhnd>0)  (73) 

J -m  u «.-*  ° 

when  the  path  of  integration  does  encircle  the  origin. 

The  terms  with  inverse  powers  of  6 are  monotonic  while  the  term  with  e~*  is 
oscillatory.  The  range  of  0 over  which  the  oscillatory  term  should  be  included  requires 
a determination  of  the  value  of  0 for  which  Jrn  6 = 0 Insofar  as  the  exponential  integral 
is  an  analytic  function  of  0.  the  path  of  integration  in  the  d-plane  may  be  deformed. 
The  point  where  the  path  of  integration  crosses  the  real  axis  may  be  displaced  to  a 
remote  position  on  the  real  axis  where  the  oscillatory  term  is  negligible.  The  oscillatory 
term  may  be  excluded  for  all  0 when  x is  positive. 


Monotonic  Terms 

The  monctonic  contribution  to  the  potential  is  given  by  the  approximation 

* —I  r Jd 


* n!  f cos*n0  d$ 

(74i 

and  the  monotonic  contributions  to  the  derivatives  of  are  given  by  the  approximations 


• n> 
i L (-0n!L  1 

»-i  n J 

r cos*"*1®  dd 

P (M  + tai)n*‘ 

(75) 

i S <-0"~ 

*%•!  n J 

f cos*""*®  sin  ®cf® 

P (m  + vw)nM 

(76) 

- L <-i)"—  d 

n-1  n J 

r cos *n-*®d® 

P ( n + iu)n *' 

(77) 

The  integrals  of  powers  oi  n + iu  are  computed  and  are  stored  in  an  array.  The  powers 
of  cos®  are  expressed  in  terms  of  powers  of  m +■  icj  The  results  are  combined  to  give 
the  terms  in  the  derivatives  of 

Let  r.  0 be  polar  coordinates  in  physical  space  The  coordinates  are  related  by  the 
equations 


x » r cos  p 


V **  r sin  0 


Then  n + iu  is  given  by  the  equation 


fi  + iu  m ft  + % %/**  + y*  co *(0  - 0)  (70) 

The  integral  of  the  reciprocal  of  fx  + iu  is  given  by  Formula  300  on  page  41  of  Peirce's 
Tabl*  of  Integrals'*.  The  basic  integrals  of  powers  of  m + is*  are  given  by  the  equations 


~ <£d5-  1 

2w  J 

(80) 

Vx*  + y*  J*  dO  _ Vx*  + y* 

2 If  J JA  4*  UJ  \Zxf  + y*  -f  j** 

(81) 

A minus  sign  is  required  for  the  second  integral  because  n is  negative.  Integrals  of 
higher  and  lower  order  can  be  generated  through  integration  and  differentiation. 

Let  the  parameter  u be  defined  by  the  equation 


u 


M 


(82) 


and  let  the  integral  /„(u)  be  defined  by  the  equation 


A.(u) 


JL  <£  (m  + <«)w 

(**  + y*)* 


dO 


(83) 


If  n is  positive,  the  integrals  are  generated  by  the  recurrence  equation 


Wu)-(n+l)  £*/,«)  dt+{ 


(~  l)"*'(n  + 1)! 


(n  even) 


(n  odd) 


(84) 


where  the  constant  of  integration  is  derived  from  the  beta  function,  and  is  nonzero 
only  when  n is  odd.  Each  integral  is  represented  by  a power  polynomial,  for  which  the 
polynomial  of  lowest  degree  is  given  by  Equation  (80)  The  polynomials  are  integrated 
term  by  term  with  the  aid  of  Equation  (3)  in  Appendix  fl.  If  n is  negative,  the  integrals 
are  generated  by  the  recurrence  equation 

/n.,(u)  ^ ~ ~ /,(u)  (85) 

rv  au 

Each  integral  is  represented  by  the  quotient  between  a power  polynomial  and  a radical, 
for  which  the  quotient  of  highest  degree  is  given  by  Equation  (81).  The  quotients  are 
differentiated  term  by  term  with  the  aid  of  Equation  (6)  in  Appendix  B. 
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Oscillatory  Term* 

The  oscillatory  contribution  to  the  potential  is  given  by  the  approximation 


^,(x.  y.  *)  — 4i  J*  s~*d< 


for  which 


1 


(100) 


(101) 


is  a monotonic  factor  of  the  integrand,  and  the  oscillatory  contributions  to  the 
derivatives  of  are  given  by  the  approximations 


for  which 


Vl  + f*  •"*  dt 

J — m 

(102) 

J fv  1 + <*  t~*  dt 

(103) 

(104) 

tVTT? 

i + f* 

(105) 

vT+T* 


are  monotonic  factors  of  the  integrands. 

Required  for  Taylor  series  expansions  of  the  monotonic  factors  are  derivatives  which 
can  be  expressed  as  the  quotients  between  power  polynomials  in  t and  powers  of 
V 1 + <*.  The  algorithm  for  generating  successive  derivatives  of  the  quotients  is  given 
by  Equation  (6)  in  Appendix  B 

The  series  converge  in  the  f-plane  only  within  that  circle  which  is  centered  at  the 
center  of  expansion  and  passes  through  the  nearest  singularity  at  ±i.  When  each  series 
of  limited  convergence  is  multiplied  by  an  exponential  function  of  its  argument  and 
the  proouct  is  integrated  from  -«  to  the  seriss  becomes  asymptotic. 

Evaluation  of  the  integrals  can  be  completed  if  the  variable  t is  replaced  by  a new 
variable  u such  that  the  parameter  £ is  a polynomial  of  at  most  the  third  degree  in 
the  variable  u 

The  parameter  S and  its  derivatives  are  given  in  terms  of  t by  the  equations 


dS 

dt 

d*6 
dt * 


- ,tt(l  + t*)  - i(x  + yt)V 1 + f* 
,(y  + r<  + 2yf*) 


- 2 m<  - v 

-2,1  - i 


Vl  + f* 

(x  + 3yf  + 2y<*) 

(1  +«*)* 
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(106) 

(107) 

(108) 


iiffttiiirifirti 
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X 
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(108) 
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The  parameter  d and  its  derivatives  are  given  in  terms  of  u by  the  equations 


d 


■ <5,  + dau  + 


d;‘u»  + do"u» 
2*  + 3' 


di 

du 


i'a  +•  ^gU  + 


d^’u* 

2! 


(no) 


d*d 

du* 


do  + dg'u 


(111) 


Equivalence  of  the  expressions  defines  a series  expansion  of  f in  terms  of  u provided 
dd/du  » 0 at  the  same  value  of  d as  dS/dt  » 0 Otherwise  df/du  would  be  infinite  where 
dd/df  * 0.  The  equation  dd/df  - 0 has  roots  where  the  values  of  f are  ft  and  f(  and 
the  values  of  d are  d,  and  dt. 

Near  tl.e  centerline  behind  the  source  the  points  of  stationary  phase  are  far  apart. 
The  path  of  integration  in  the  f-plane  may  be  deformed  to  pass  through  the  points 
f,  and  f,  of  stationary  phase.  The  contribution  to  integration  is  large  near  t,  and  ft. 
and  is  small  between  f,  and  ft.  Each  integration  from  f,  or  ft  to  an  intermediate  point 
may  be  approximated  by  an  integration  to 

The  polynomial  for  d is  simplified  by  setting  dj,”  = 0 The  conventional  substitution 

u-u-^7,  (112) 

eliminates  also  di.  The  constant  da  is  given  by  the  uubstitution 

d#-d  (113) 


and  the  second  coefficient  6 i is  given  by  tne  substitution 


(114) 


r 


i 


i 

\ 

\ 


i 


for  which  f has  the  values  f,  or  f,  where  dd/df  = 0.  The  algorithm  for  converting  the 
variable  of  integration  from  f to  u is  given  by  Equation  (15)  in  Appendix  B. 

Integration  of  * with  respect  to  u is  given  in  terms  of  exponential  functions  by 
the  equation 


(115) 


Integration  of  the  product  of  u and  e~ * is  zero  by  symmetry  Integration  of  the  product 
of  a power  of  u and  a~*  is  given  by  the  recurrence 


utn  *e~*  du 


or  after  iteration  by  the  equation 


du 


(2:.)?(3w)» 

2-n!(di')n+* 


(116) 


(117) 
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Insofar  as  integration  gives  monotonic  terms,  summation  can  be  continued  until  a 
new  terra  is  larger  than  the  preceding  terra. 

Along  the  line  of  wave  crests  the  points  of  stationary  phase  are  close  together  The 
path  of  integration  in  the  (-plane  may  be  deformed  to  pass  midway  between  the 
points  (,  and  f(  of  stationary  phase. 

The  quadratic  term  in  the  cubic  polynomial  can  be  eliminated  by  the  conventional 
transformation 


i*  •*  u - 


«s 

*e" 


<U8) 


The  variable  u thus  can  be  redefined  so  that  dj(  ■ 0.  Inasmuch  as  u is  a new  variable, 
it  may  be  scaled  arbitrarily.  Let  it  be  scaled  so  that  dj'  • -2f.  If  the  limits  of  integration 
lie  above  the  complex  value 


l*«slt 

N 


(119) 


then  the  real  part  of  -d  goes  to  in  the  l*mit  as  N - The  integral  of  e“4  with 
respect  to  u is  guaranteed  to  converge  if  the  path  of  integration  lies  just  above  the 
real  axis 

The  equation  dS/du  « 0 is  . mplified  to  the  equation 


u*--idi  (120) 

for  which  the  value  of  d is  given  by  the  expression 

d,±J«di)»  (121) 

The  values  of  da  and  d#  are  given  in  terms  of  the  values  of  <5,  and  dt  by  the  equations 

d.  - $<d,  + dt)  (122) 

di«-il!(di-d,)}*  (123) 

The  value  of  ( for  which  d = d0  is  determined  by  Newton-Rephson  iteration.  The  value 
of  ( at  which  the  iteration  is  started  is  given  by  the  approximation  ( ~ *((,  + (t).  The 
algorithm  for  converting  the  variable  integration  from  ( to  u is  given  by  Equation 
(17)  in  Appendix  B 

Integration  of  e*4  with  respect  to  u is  given  in  terms  of  Airy  functions  by  the 
equation 

J*  e'4  du  * 27T«"4°At(idi)  (124) 

Integration  of  the  product  of  a power  of  u and  *~4  is  given  by  the  equation 

J une-4du- 2rrt‘*8(-i)"At<n>(idi)  (125) 

where  Ai<n)(id|,)  is  the  nth  derivative  of  Ai(tdi)  The  algorithm  for  generating  the 
successive  derivatives  is  given  by  Equations  (26)  and  (27)  in  Appendix  B. 
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Inasmuch  as  integration  gives  terms  which  alternate  between  two  series  of  values, 
summation  is  continued  until  a new  term  is  larger  than  both  of  the  two  previous 
terms  or  until  the  terms  make  no  change  in  the  sum.  In  a range  of  arguments  the 
terms  increase  at  first  before  they  decrease,  and  finally  increase.  In  order  to  keep 
the  initial  increase  from  triggering  a premature  termination  of  the  series,  the  first 
three  terms  are  evaluated  prior  to  the  start  of  sensing  for  termination  of  the  series. 


Accuracy  and  Efficiency 

Transition  from  the  quadratic  approximation  to  the  cubic  approximation  is  placed 
where  the  errors  of  computation  are  the  same  for  both  approximations.  The  cubic 
approximation  is  used  where  the  coordinates  satisfy  the  inequality 


M > *!vl 


(126) 


Transition  from  truncation  error  to  rounding  error  occurs  at  a critical  radius  from 
the  source.  The  limiting  radius  for  acceptable  truncation  error  varies  with  direction. 
The  critical  radius  is  given  by  the  empirical  equation 


+ ^ 2 7 + 21 


X*  4-  y* 


(127) 


The  rounding  error  for  radii  larger  than  the  limiting  radius  also  varies  with  direction. 
The  rounding  error  along  the  critical  line  of  wave  crests  increases  with  decrease  in 
the  depth  below  the  free  surface.  The  relatively  large  errors  on  the  critical  line  are 
attributed  to  the  close  approach  to  each  other  of  the  points  of  stationary  phase.  Thus 
a perturbation  of  the  position  of  the  points  of  stationary  phase  causes  a significant 
change  in  the  error  in  velocity.  The  error  is  tolerated  insofar  as  the  asymptotic 
approximation  is  faster  than  the  integration  by  parts.  Errors  and  times  are  illustrated 
by  dotted  curves  in  Figures  5 and  6 


INTEGRATION  BY  PARTS 

Interpolation 

In  the  integration  with  respect  to  t the  integrands  are  the  products  of  the  monotonic 
factor 


and  the  oscillatory  factor 


T - 

J—  u 


du 


(128) 


(129) 


or  the  integrands  are  the  products  of  the  monotomc  factors 

Vi  + 1*  fW  + t* 

and  the  oscillatory  factor 


'J1S 


du 


1 + t*  (130) 

(131) 
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If  the  monotonic  factors  and  the  oscillatory  factors  are  expressed  in  terms  of  a 
common  argument,  then  the  integration  can  be  completed  through  integrations  by 
parts.  The  parameter  5 is  expressed  by  the  equation 

i-rt  + t (132) 


where  i)  is  a center  of  expansion  and  t is  the  common  argument  of  the  monotonic 
factors  and  the  oscillatory  factors. 

The  real  part  of  rj  always  is  positive.  The  real  part  of  tj,/*  always  is  negative,  because 
rj  is  on  the  same  branch  of  the  Riemann  surface  as  -«■.  The  sign  of  ct/(  is  arbitrary. 
The  sign  is  adjusted  so  as  to  make  the  phase  angle  of  t1/ * equal  to  half  the  phase 
angle  of  c. 

Subtraction  of  rj  leads  to  the  equation 

* “ - (h  + M)  ~ fit*  - i(x  + yf)V  1 + t*  (133) 

* 

and  division  by  f*  leads  to  the  equation 


Taylor  series  expansions  express  c directly  in  terms  of  t,  then  transformations  of 
variable  express  t conversely  in  terms  of  t.  The  algorithms  for  the  series  expansions 
and  the  transformations  of  variable  are  the  basis  for  a previous  subroutine  which  has 
been  described  in  a previous  report.  Tnis  subroutine  spent  too  much  time  in  processing 
the  Taylor  series  expansions.  A gain  in  speed  by  a factor  of  seven  has  been  achieved 
when  the  Taylor  series  expansions  have  been  replaced  by  Lagrange  interpolation.  The 
algorithms  for  the  Lagrange  interpolation  are  the  basis  for  the  current  subroutine 
which  is  described  in  the  present  report. 

Series  expansions  of  the  monotonic  factors  are  derived  by  Lagrange  interpolation 
between  eleven  discrete  values  of  the  monotonic  factors.  The  discrete  values  are 
computed  at  each  discrete  argument  in  a set  whose  spacing  between  arguments  is 
proportional  to  the  spacing  between  the  roots  of  a Chebyshev  polynomial.  The  accuracy 
of  approximation  between  discrete  arguments  then  tends  to  be  uniform  in  accordance 
with  the  Chebyshev  criterion.  Where  the  path  of  integration  is  curved,  the  spacing 
along  the  curve  is  approximated  by  chords  which  span  segments  of  the  curve.  The 
range  of  approximation  is  determined  by  the  configuration  in  the  (5-plane,  while  the 
path  of  integration  is  specified  hi  the  f-plane.  The  position  in  the  f-plene  for  a 
specified  chord  in  the  4-plane  is  derived  by  Newton- Raphson  iteration. 

The  natural  path  of  integration  is  the  real  axis  in  the  f-plane.  but  this  path  passes 
between  two  singularities  where  d<5/cff  = 0.  Many  intervals  of  integration  would  be 
required  in  the  vicinity  of  the  singularities.  These  can  be  reduced  to  a single  interval 
near  the  first  singularity,  and  they  can  be  reduced  to  fewer  intervals  near  the  second 
singularity,  if  the  path  of  integration  is  displaced  to  pass  right  through  the  first 
singularity.  The  integration  is  completed  in  a single  final  interval. 
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In  the  integrations  with  respect  to  t the  monotonic 

factors  are 

1 

di 

(135) 

dt 

and 

n/1  + t * 

tVi  + 11 

(1  -M  *) 

(136) 

di 

di 

di 

dt 

dt 

dt 

Different  expansions  are  required  for  each  interval  of  integration. 

In  the  first  interval 

of  integration,  the  center 

of  expansion  r)  is  located 

at  that 

value  of  i where  t ■ ft, 

and  di/dt  = 0.  The  parameters 

1 

t* 

di 

(137) 

dt 

and 

1 

1 

I 

tWl  + t 2 

**<vi  + 1* 

C*(l  +f*) 

(138) 

di 

di 

di 

dt 

dt 

dt 

are  computed  for  a discrete  jet  of  values  of  e,/* 

Their  limiting  values  at  i = 

0 are 

1 

1 

(139) 

and 

VI  + t,* 

f.vi  + «,* 

(1  ♦ *,*) 

(140) 

i 

i 

1 

\ dt * / 

\dt*J 

\IU‘I 

Expansion  in  series  leads  to  the  representation  of  the  monotonic  factors  as  ascending 

series  in  powers  of  c1  •* 

which  begin  with  € 1 2 

In  the  continuation  of  integration,  the  center  of  expansion  tj  is  located  at  the  center 
of  the  range  of  approximation.  The  monotonic  factors  are  expressed  directly  as  series 

in  powers  of  e. 

In  the  final  interval 

of  integration,  the  center 

of  e 

rpansion  tj  is  located 

at  that 

value  of  S where  di/dt 

= 0 and  5 is  nearest  to  + <*>  The  parameters 

i 

e* 

i di 

(141) 

dt 

and 

Vi  + f* 

fVl  + f* 

(1  + f*) 

(142) 
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i di 

dt 

* dt 

£ dt 
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are  computed  for  a discrete  set  oi  values  o i Their  limiting  values  at  e - <•>  are 


(-/x?iy)* 


(“  M * iy)2 


(-/XT  iy)* 


where  the  ± signs  are  determined  according  to  whether  the  direction  of  integration 
is  ± in  the  t- plane.  Expansion  in  series  leads  to  the  representation  o'  the  monotonic 
factors  as  descending  series  in  powers  of  £~l/i  which  begin  with  c*1'1*. 

The  range  of  expansion  of  the  Lagrange  interpolation  is  subject  to  limitations  which 
are  similar  to  the  limitations  on  the  range  of  expansion  of  the  Taylor  series.  The 
expansions  are  limited  to  the  primary  branch  in  the  6 -plane  as  long  as  the  radius 
of  expansion  is  less  than  the  distance  from  the  center  of  expansion  to  the  points 
where  <3  = 0 and  t = ± i.  The  expansions  are  conv  ent  if  the  radius  of  expansion  also 
is  less  than  the  distance  from  the  ceui.fi  oi  expansion  to  the  nearest  point  in  the 
primary  branch  where  dd/dt^O  The  actual  range  of  approximation  is  a fraction  of 
the  radius  of  convergence.  The  value  of  the  fraction  is  determined  by  computation. 
In  the  initial  interval  of  integration  the  fraction  for  expansion  in  terms  of  €1/z  is  one 
half,  in  the  continuation  of  integration  the  fraction  for  expansion  in  terms  of  c is  two 
thirds,  and  in  the  final  interval  of  integration  the  fraction  for  expansion  in  terms  of 
£~1/l  is  one  half.  Boundaries  of  approximation  and  convergence  are  illustrated  by 
dotted  lines  in  Figures  3 and  4 for  the  initial  interval  of  integration. 

During  the  stepping  from  one  interval  to  the  next  interval  the  center  of  expansion 
for  the  new  interval  is  estimated  on  the  assumption  that  the  new  interval  is  a linear 
extension  of  the  old  interval.  The  position  of  the  center  of  the  new  interval  is  estimated 
by  an  application  of  the  cosine  law  to  a triangle  with  vertices  at  the  old  center  of 
expansion,  the  nearest  singularity,  and  the  new  center  of  expansion 


Recurrence 

The  terms  of  the  ascending  series  are  integrated  with  the  aid  ci  the  recurrence 
equation 

r e1  r*  un  rn+*  e‘ 

une -(,+•)  — dfdu*  du-  «»« -<*+*>  — dt 

Jo  J — i JoV+u  t 

+ nJun-ie'(^u)£U~dtdu  (145) 

which  is  started  with  the  aid  of  the  special  equations 

1 "x  JI  t * = 2e'”  L 3 un'iv dt  (u6) 

rv<->  r^u-iogfi^w*  r ru  w 


and  with  the  aid  of  tha  racurrenca  equation 
r>  r+* • r*  „■*- 


I u-e  £*  ?;  « du  - n £ £1  du  - iV-»  £*  ^ df 

+ n £ u»->*-<e*w)  J'  " ^ df  du 


C?--G 

r- 

(149) 

r***  -i 

L 7* 

(150) 

which  ia  started  with  the  aid  of  the  special  equations 

<**,-<»•*«*>  t«  2 /€>l 


J*  * ^ ' ji  d*  du  ~ f Jdt  ~ * Tdt  (150) 

The  recurrence  is  cycled  in  the  reverse  direction  when  |«|  S 7, 

Required  for  ascending  senes  are  integrals  which  are  generated  by  the  recurrence 
equation 

r vn  e*  f*  u*'1 

I du  = tj  I du  (151) 

Jo  V + u n Jo  »)  + u 

The  recurrence  is  started  with  the  special  integrals 

f -A -■%•*(;)*  052) 

*'•  u*(q  + u)  q*  ',T" 


r,i  \n/ 

(152) 

10,(’%) 

(153) 

The  integrals  are  given  by  the  series  expansion 


Jo  + ■**  U 77  *-on  + * + 1^/ 


(Id  < Ini)  054) 


The  recurrence  (151)  removes  progressively  the  terms  of  lowest  order  from  the 
expansion,  whence  the  residue  does  not  retain  the  terms  of  highest  order  with  enough 
accuracy  when  !e|  <<  |rj|. 

The  expansion  (154)  is  generated  in  the  direction  of  descending  order  by  the 
recurrence 


rjfndu, ].[!:_  r_MLdj 

Jo  o + « it"  Jo  n + u J 
which  may  be  started  with  the  limiting  approximation 

r u*  e**1 

Jo  V + u (,V  + I X77  + e) 


(,V  - »)  (156) 


In  view  of  the  convergence  condition  for  the  series  expansion,  the  recurrence  (155) 
can  be  used  only  when  ie|  < !*?(. 
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The  terms  of  the  descending  series  are  Integrated  with  the  aid  of  the  recurrence 
equation 


r r 4 * 

J,  J — < njt  q+u  m*  J_.  t 

-ij  «-*.—>  J yrftctu 

which  is  started  with  the  aid  of  the  special  equations 

r ^ n-~ 

and  with  the  aid  of  the  recurrence  equation 

r— r;  ? - 

-if-— £' 

which  is  started  with  the  aid  of  the  special  equations 

£ uKl***>  £ dt  ctu  - £ 'j  dt  + *-’£  ^tan-'^^df  (161) 

rr>*~  a1  n*’  e1 

L ?**- -^’L  idt 

r ^ r.  ?d<du="^°*(i+7) + -d*h)* 


(157) 

(156) 

(159) 


(160) 


(162) 

(163) 

(164) 


The  recurrence  is  cycled  in  the  reverse  direction  when  ;ef  2 27. 

Required  for  descending  series  are  integrals  which  are  generated  by  the  recurrence 
equation 

j ■**"* 


r^!du=i[_L.rjLidj 

J,  n + u nine"  J.  v + u 


The  recurrence  is  started  with  the  special  integrals 

du  2 

i “ 

*(v  + u)  r) 


^ U*(q  + u)  »)*  ' ' 

faSw-i-H) 


(165) 

(166) 
(167) 
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(W  < 1*1)  (i««) 


j 


i 
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The  integrals  are  given  by  the  eeriee  expansion 


u~* 
f)  + u 


du-^E 

" »•« 


(-0* 
n + k 


The  recurrence  (165)  removes  progressively  the  terms  ol  highest  order  from  the 
expansion,  whence  the  residue  does  not  retain  the  terms  of  lowest  order  with  enough 
accuracy  when  !»)|  <<  |«| 

The  expansion  (168)  is  generated  in  the  direction  of  ascending  order  by  the  recurrence 


du  “ 


1 

tu* 


tj  + u 


dti 


(169) 


which  may  be  started  with  the  limiting  approximation 


r 


ij + « 


du  ~ 


(N-  l)(f|  + «) 


(Af -»  -)  (170) 


In  view  of  the  convergence  condition  for  the  series  expansion,  the  recurrence  (169) 
can  be  used  only  when  |tjl  < |*|. 

The  complex  exponential  integral 


(171) 


and  the  complex  Fresnel  integral 


(172) 


are  obtained  from  subroutines*4.  It  is  necessary  to  apply  corrections  to  bring  the  path 
of  integration  below  each  singularity  at  t » 0.  The  corrections  are  applied  when 


Ant;  > 0 

For  the  exponential  integral  the  correction  is 

+ Ziri  e~* 

but  for  the  Fresnel  integral  the  correction  is 

- zVnit~*  - 2e~* 


j—  it 


The  logarithmic  exponential  integral 


(173) 

(174) 

(175) 

(176) 


I 


| 


(j 

i 

i 

| 
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and  the  arctangential  Fresnel  integral 


Ulna, jfr* 


are  obtained  from  subroutines**  H i*  neceaeary  to  apply  correction*  to  bring  the  path 
of  integration  below  each  singularity  at  ( » For  the  logarithmic  exponential  integral 
the  correction  it 


and  for  the  arctangential  Fresnel  integral  the  correction  is 


In  the  evaluation  of  the  ascending  series  the  arctangential  Fresnel  integral  is  corrected 
only  when 


jUti  0 and  frn « > 0 or  Jt*  * < 0 and  fn 


!n  the  evaluation  of  the  descending  series  the  integrals  ore  corrected  only  when 

Jtec>0  and  /mc>0  (181) 

When  the  differences  are  taken  between  values  for  the  upper  and  lower  limits  of  t the 
correction  for  the  logarithmic  exponential  integral  is  reduced  to 


- 2rrt* 


and  the  correction  for  the  arctangential  Fresnel  integral  is  reduced  to 

r*  *' 

- we*’1  -r  dt  (183) 

J—  *1 

The  arctangential  Fresnel  integral  and  the  arccotangential  Fresnel  integral  are  related 
in  accordance  with  the  equation 

•’£  ft  “""(j)1  “ - 1 •-£  fi " - •-£  'j  ‘“"(i)' " (,84) 

which  is  used  in  connection  with  the  evaluation  of  the  ascending  series 


Accuracy  and  Efficiency 

The  integration  by  parts  gives  acceptable  accuracy  and  efficiency  everywhere.  Errors 
and  times  are  illustrated  by  solid  curves  in  Figures  5 and  6. 


PROGRAMMING 
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Trapezoidal  integration  for  velocity  pnlential  ia  provided  by  Subroutine!  CKPSVP  and 
CEXPU.  while  aaymptotic  approximation  and  Integration  by  parts  are  provided  by 

SUBROUTINE  pswtvp  (ak,  ah.  ax.  ay.  az,  FP)  j 

FORTRAN  SUBROUTINE  FOR  VELOCITY  POTENTIAL  OF  POINT  SOURCE  ( 


The  critical  wave  number  ia  given  in  argument  AK.  and  the  depth  h of  the  source  i 

ia  given  in  argument  AH.  The  coordinates  x.  y,  z of  a point  are  given  in  arguments  j 

AX.  AY,  AZ.  The  velocity  potential  if  is  placed  in  function  FP.  Calls  are  made  to  auxiliary 
subroutines  which  are  listed  in  Table  I.  j 

Trapezoidal  integration  for  velocity  field  is  provided  by  Subroutines  CKPSVP  and  j 

CEXPDf.  while  asymptotic  approximation  and  integration  by  parts  are  provided  by 

SUBROUTINE  PS'VTVi  (AK.  AH.  AX.  AY,  AZ,  FU.  FV.  FW)  j 


FORTRAN  SUBROUTINE  FOR  VELOCITY  FIELD  OF  POINT  SOURCE 


The  critical  wave  number  is  given  in  argument  AK,  and  the  depth  h of  the  source 
is  given  in  argument  ah.  The  coordinates  x.  y.  x of  a point  ore  given  in  arguments 
AX,  AY,  AZ.  The  components  u,  v,  vu  ct  velocity  are  placed  in  functions  FU,  FV.  FW.  Calls 
are  made  to  auxiliary  subroutines  which  are  listed  in  Table  I. 


TABLE  1 

AUXILIARY  SUBROUTINES 

Name  Function 

C7CTRM Performs  complex  vector -matrix  multiplication. 

CLCRNP Synthesizes  complex  Lagrangian  polynomials. 

C8SSLK Computes  the  modified  Bessel  function  of  the  second  kind. 

CLGMCI  Computes  the  complex  logarithmic  integral. 

CEXPLI Computes  the  complex  exponential  integral. 

CFRNII Computes  the  complex  Fresnel  integral. 

CEXEXI Computes  the  complex  exponential  exponential  integral. 

CEXFRI Computes  the  complex  exponential  Fresnel  integral. 

CLGEXI Computes  the  logarithmic  exponential  integral 

CATFRI Computes  the  arctangential  Fresnel  integral. 
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DISCUSSION 


The  fundamental  formulation  of  the  Green  function  can  be  transformed  into  other 
formulations  in  which  the  added  Fourier  integrals  have  various  amplitudes.  Methods 
of  evaluation  by  contour  integration  in  the  complex  plane  of  a coordinate  in  wave 
number  space  have  been  published  in  the  literature.  The  final  evaluation  in  each 
method  is  an  application  of  a quadrature  rule  to  the  azimuthal  integration.  Any 
systematic  quadrature  rule  other  than  the  trapezoidal  rule  for  equally  spaced  arguments 
is  equivalent  to  a combination  of  trapezoidal  rules  of  reduced  order.  It  can  have 
superior  accuracy  only  by  accident. 

Radial  integration  along  a complex  path  was  suggested  originally  by  Pond’  in  1957. 
and  was  programmed  for  computation  by  DiDonato7  in  1958.  In  the  DiDonato  method, 
the  amplitude  had  its  original  form  as  expressed  by  the  equat'on 


A(k.  8)  = 


1 (x0  + X COS*g) 

2jtx  (x0  - x cos*0) 


(165) 


Radial  integration  was  in  the  complex  plane  of  the  radial  coordinate.  Contour  integration 
and  evaluation  of  the  residue  of  a pole  led  to  the  sum  of  a single  integral  and  a double 
integral.  The  single  integral  had  a variable  upper  limit.  The  radial  integration  was  by 
Simpson  rule,  and  the  azimuthal  integration  was  by  Gauss  rule. 

Radial  integration  along  a complex  path  was  suggested  again  by  Smith.  Giesing.  and 
Hess*  in  1963,  and  by  Webster*  in  1969.  It  was  programmed  for  computation  by  Adee* 
in  1973.  In  the  Adee  method,  the  amplitude  is  resolved  as  expressed  by  the  equation 


A(k.  8) 


1 

2 n* 


+ 


*o 

trx(x8  - K COS*0) 


(186) 


The  flrst  term  of  the  expansion  is  the  amplitude  of  a negative  image  source.  Radial 
integration  in  wave  number  space  was  performed  along  a complex  path  where  the 
integrand  is  not  oscillatory.  Radial  integration  was  by  series  expansion  and  Simpson 
rule,  then  azimuthal  integration  was  completed  by  Simpson  rule.  The  radial  integration 
evaluated  just  the  conventional  exponential  integral,  and  the  azimuthal  integration 
was  applied  to  the  same  integrand  as  in  the  case  of  the  trapezoidal  rule. 

A formulation  of  the  Green  function  has  been  given  by  Gadd10  in  1970.  In  the  Gadd 
method,  the  amplitude  is  resolved  in  accordance  with  the  alternative  equation 


A(k.  8) 


1 cos*0 

2jt x it (k9  -k  cos *0) 


(187) 


The  flrst  term  of  the  expansion  is  the  amplitude  of  a positive  image  source.  As  k -*  °° 
in  this  method  the  integrand  does  not  diminish  as  fast  as  in  the  conventional  method. 
If  the  source  is  close  to  the  surface,  a large  part  of  the  integration  is  devoted  to  the 
computation  of  twice  the  potential  of  the  negative  image.  Slow  convergence  of  the 
integrand  makes  accurate  evaluation  difficult. 

That  radial  integration  evaluates  an  exponential  integral  was  recognized*1**  in  1959 
The  exponential  integral  was  used  in  subroutines  which  were  reported16  in  1965.  That 
direct  evaluation  of  the  exponential  integral  is  more  efficient  than  the  conventional 
radial  integration  was  appreciated  by  Shen  and  Farell11  in  1975.  In  the  Shen  and  Farell 
method,  the  exponential  integral  is  evaluated  by  reference  to  a subroutine.  In  the 
azimuthal  integration,  the  angle  8 is  replaced  by  its  tangent  t.  The  azimuthal  integration 
is  completed  by  Simpson  rule.  The  computer  is  required  to  halve  the  interval  of 
integration  until  the  integral  remains  constant  to  within  a tolerance.  Except  for  a 
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factor  d8/dt,  the  integrand  is  the  same  as  the  integrand  for  the  trapezoidal  integration. 

The  integration  in  wave  number  space  has  been  analysed  in  terms  of  Cartesian 
coordinates  by  Anderson12,15  and  by  Noblesse14  in  1976.  In  the  Noblesse  method, 
contour  integration  in  the  complex  plane  of  the  longitudinal  Cartesian  coordinate 
leads  to  a resolution  of  the  Green  function  into  a monotonic  term,  which  is  symmetric 
fore  and  aft,  and  an  oscillatory  term,  which  trails  behind  the  source.  The  monotomc 
term  has  an  exponential  integral  in  its  integrand,  but  the  argument  of  the  exponential 
integral  is  not  the  same  as  the  conventional  argument.  The  real  part  of  the  argument 
changes  sign  but  the  imaginary  part  remains  negative  during  azimuthal  integration. 
The  argument  of  the  exponential  integral  goes  to  zero  at  the  limits  of  integration, 
and  the  exponential  integral  presents  there  a pair  of  logarithmic  singularities. 
Subtraction  of  a logarithm  from  the  integrand  gives  a function  more  suitable  for 
quadrature,  then  the  logarithm  is  restored  with  analytic  integration.  The  oscillatory 
term  is  just  twice  the  conventional  single  integral,  and  presents  the  same  problem  of 
evaluation.  Although  the  sum  of  the  monolonic  term  and  the  oscillatory  term  is 
continuous,  the  individual  terms  have  discontinuities  in  derivative  at  the  transverse 
plane  through  the  source.  They  do  not  satisfy  the  free-surface  boundary  conditions 
at  this  transverse  plane. 

An  attempt  to  compare  the  present  method  with  the  various  other  methods  has 
been  thwarted  by  an  inability  to  make  the  programming  of  the  other  methods  run 
properly  on  the  computer  in  this  laboratory. 

An  effort  has  been  made  to  interpolate1*  in  a table  of  velocities.  A table  of  36900 
entries  was  computed  on  the  Naval  Ordnance  Research  Calculator  when  that  computer 
was  waiting  to  be  demolished.  The  interpolation  was  founded  on  the  assumption  that 
the  velocities  could  be  expressed  as  the  product  of  a monotonic  function  and  the 
exponential  function  of  a monotonic  argument  This  is  true  where  the  asymptotic 
approximation  is  valid,  but  close  to  the  source  the  divergent  wave  system  and  the 
transverse  wave  system  are  inextricably  interlocked.  Efforts  to  interpolate  had  only 
indifferent  success.  In  order  to  interpolate  the  table  would  have  to  be  so  large  as  to 
be  unwieldy.  The  three-way  interpolation  could  be  as  costly  as  the  direct  evaluation. 

Measurements  of  wave  height  of  a Rankine  ovoid”"27  confirm  the  validity  of 
computations  of  wave  height  by  subroutine.  The  Rankine  ovoid  is  that  streamline  which 
is  generated  by  a simple  source  and  sink  on  a line  parallel  to  the  free  stream  A model 
of  the  Rankine  ovoid  has  been  towed  in  the  model  basin  and  wave  heights  have  been 
measured  by  Shaffer25  The  measured  wave  heights  agreed  in  phase  with  the  computed 
wave  heights  but  were  smaller  in  amplitude.  The  model  was  towed  by  a long  stranded 
cable,  which  undoubtedly  set  up  a fully  turbulent  boundary  layer  on  the  model. 
Inasmuch  as  vorticity  in  the  boundary  layer  would  meet  partially  the  boundary 
conditions  on  the  boundary  of  the  ovoid,  it  was  necessary  for  the  source  distribution 
to  meet  only  partially  the  boundary  conditions.  The  agreement  between  measured  and 
computed  wave  heignts  is  excellent  in  view  of  the  experimental  difficulties. 


CONCLUSION 

The  trapezoidal  integration  is  efficient  in  the  near  field  of  a deep  source,  but  is 
inefficient  in  the  far  field  of  a shallow  source.  The  asymptotic  approximation  is  efficient 
far  from  the  source,  but  is  inaccurate  along  the  critical  line  of  wave  crests.  The 
integration  by  parts  is  most  uniformly  accurate  and  efficient  at  all  depths 
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TRANSIENT  WAVE 


\ 
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tn  the  analysis  of  a transient  wave,  the  uniform  motion  of  a source  along  a horizontal 
line  is  simulated  by  a succession  of  pulses  at  equal  intervals  of  distance  and  time. 

Let  *.  v.  z,  f be  coordinates  and  time  with  origin  of  coordinates  at  the  surface  of 
the  fluid  initially  at  rest.  Let  a fixed  source  be  created  at  a distance  ti  below  the 
origin.  The  velocity  is  the  negative  gradient  - of  a velocity  potential  9,  which  is  a 
solution  of  Lap'ace’s  equation.  The  dynamical  equation  is  the  complete  Bernoulli 
equation. 


+ iOV)*  + ~ gz  » constant 

dt  P 


(1) 


where  p is  the  density,  p is  the  pressure,  and  g is  the  acceleration  of  gravity.  Let  the 
surface  of  the  fluid  be  given  by  the  equation 


* + <(*.y.  0 = o 

where  ( is  the  elevation  of  the  surface.  The  linearized  Bernoulli  equation  is 

0a 

- ZZ.  + g{-  = 0 

at 

at  the  free  surface,  and  the  kinematic  equation  is 

a<  n 

+ — * 0 

oz  at 

Elimination  of  ( from  Equations  (3)  and  (4)  leads  to  the  boundary  equation 

dt*  9 dz 


(3) 


(3) 


(4) 


(5) 


which  must  be  satisfied  by  the  transient  potential. 

At  the  instant  of  creation  of  the  source  the  fluid  is  accelerated  into  motion.  All 
terms  in  the  Bernoulli  equation  remain  bounded  at  the  free  surface.  I?  each  term  is 
integrated  with  respect  to  time,  then  the  term  d<p/dt  becomes  the  change  of  potential, 
but  the  other  terms  vanish  in  the  limit  of  instantaneous  formation  of  the  source, 
the  potential  remains  constant,  and  the  boundary  conditions  require  that  the  formation 
of  a source  beneath  the  surface  be  accompanied  by  the  formation  of  an  image  over 
the  surface. 

The  potential  of  the  surface  disturbance  is  given  by  the  equation 
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<p(x,  y.  z.  f)  = <px(x.  y.  z)  4 <pt(x,  y.  z)  4 <p3{x,  y,  z,  f)  (6)  l 

i 

where  is  the  potential  of  the  source  in  an  unbounded  fluid,  <pt  is  the  potential  of 

the  image  source  over  the  free  surface,  ar.d  is  the  potential  of  a transient  disturbance.  > 
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For  a unit  source  below  the  surface  the  potential  ?t  is  given  by  the  equation 
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(7) 


v**  + y*  + (*  - h)* 

and  for  the  image  source  over  the  surface  the  potential  <ft  is  given  by  the  equation 
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s»« 


Vx*  + v*  + (*  + h)* 

The  Fourier  transforms  of  the  three  potentials  are  given  by  the  equations 
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Substitution  of  potentials  in  the  free-boundary  equation  shows  that  the  Fourier 
amplitude  satisfies  the  equation 


t'A 

dt* 


+ p«4 


of  which  the  appropriate  solution  is 


A(k.  ».<)=■  — 
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If  an  instantaneous  creation  of  a finite  source  is  followed  after  a time  interval  dt  by 
the  instantaneous  annihilation  of  the  finite  source,  then  the  Fourier  amplitude  per 
unit  pulse  is  merely 


A(k.  9,  t) 


(14) 


The  potentials  and  <pt  are  zero  for  the  unit  pulse,  while  the  potential  <p3  is  given 
by  the  equation 


i rw  rm 

<f3  = - J J vpc  sin  v^f  *-«->-<«— > fa  m (15) 

Substitution  of  the  potential  in  the  boundary  equation  shows  that  the  surface  elevation 
initially  is  a hump,  but  breaks  ultimately  into  the  well  known  system  of  concentric 

waves. 

If  a unit  source  is  created  below  the  origin  at  depth  h and  then  moves  along  the 
x-axis  at  constant  speed  U,  the  potentials  and  <pz  are  restored  but  with  their 
origins  centered  over  the  source.  The  potential  <p3  me.y  be  obtained  through  integration 

2 


I 

.1 


I 

I 

I 

i 

! 


I 

I 

* 

> 

f 

t 

i 

* 


i 


i 


r 


of  the  contributions  of  pulses  which  have  been  created  during  a succession  of  differentia! 
time  intervals  dr  from  t * 0 to  r » t.  The  result  at  time  t is  given  by  the  equation 


i JJ  J Vfft  sin  y/gi{l  - T)  dKdadr 

After  integration  the  potential  is  given  by  the  equation 
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(17) 


In  the  limit  «><■*•»,  the  first  integral  vanishes  unless  cos  8 is  negative,  and  the  second 
integral  vanishes  unless  cos  8 is  positive  In  either  case  there  is  no  contribution  to 
the  integration  except  where  the  denominators  vanish.  The  substitution 


x = 
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U*co3*8 


+ ti 


(18) 


and  use  of  the  equation 
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sin  au 


du 


n 
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(<*>0)  (19) 


for  which  a is  given  by  the  equation 

a = ji/tloos  8 (20) 

converts  the  sum  of  the  first  two  integrals  into  the  single  integral  in  Equation  (34) 
of  the  text  In  the  second  two  integrals  the  integrands  are  multiplied  by  the  square 
of  a sine  whicn  is  always  positive  and  in  the  limit  as  t ■*  <■»  has  the  effect  of  weighting 
the  Cauchy  principal  value  by  | The  sum  of  the  secono  two  integrals  becomes  the 
double  integral  in  Equation  134)  of  the  text  Equation  i\~!)  is  very  instructive  insofar 
as  it  demonstrates  the  true  nature  of  the  Cauchy  principal  value 
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ALGORITHMS 


Integration  of  Sorico 

Lot  * function  4(f)  bo  expressed  by  the  equation 

-*(0  “ £ «*<* 

*-e 

The  nth  integral  4<"*(f)  io  given  by  the  definition 

4<**(1)  - Z «i"»f* 

M 

for  which  the  coefficients  are  derived  from  the  recurrence 


The  coefficients  are  stored  in  an  array  and  ere  ahifted  during  each  cycle  of  recurrence. 
The  constant  of  integration  a^)  is  replaced  by  a new  constant  in  each  cycle. 


Differentiation  of  Quotient 

Let  a function  Q(f)  be  expressed  by  the  equation 


The  nth  derivative  Q(Ki{t)  is  given  by  the  definition 


Qln)(0 


a + f*r* 


for  which  the  coefflcienta  are  derived  from  the  recurrence 

- (k  - 2n)aiV»  + (*  + 1 )«fc°  (8) 

The  coefficients  ore  stored  in  an  array  and  ere  replaced  during  each  cycle  of  recurrence. 


Tranaformation  of  Variable 

If  a dependent  variable  is  expressed  in  terms  of  either  of  two  independent  variables 
by  power  polynomials,  then  the  equality  of  the  power  polynomials  defines  implicitly 
one  independent  variable  in  terms  of  the  other,  provided  the  power  polynomials  meet 
necessary  restrictions.  Let  P be  the  dependent  variable  and  let  f,  u be  the  independent 
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variable*.  Let  the  variable  P be  given  either  by  the  equation 


or  by  the  equation 


P-  I <W 


P-  I 6,u" 


Equivalence  of  the  two  polynomial*  implies  that  t and  u are  related  by  the  equation 

t • £ e*«*  (9) 

All  of  these  polynomials  are  restricted  to  have  zero  constant  terms. 

The  mth  power  of  t is  given  by  the  definition 


tm~  £ cir’w*  (10) 

The  coefficients  cjj"*  all  are  zero  for  n < m and  the  first  nonzero  coefficient  is  given 
by  the  equation 

cjrMc.r  (a) 

The  coefficients  are  generated  by  the  recurrence 

(12) 

4-1 

The  expression  for  c^"*  contains  no  coefficient  c*l>  of  higher  order  than  , . 

The  coefficients  are  solutions  of  the  equations 

E amc<T)  - 6„  (13) 

which  express  the  equivalence  of  the  polynomials  for  P Solution  of  the  equations  is 
achieved  by  an  iteration  which  starts  with  all  coefficients  cn  set  equal  to  zero  except 
c,.  When  a,  and  6,  both  are  zero,  the  coefficient  c,  is  given  by  the  equation 

Ci  = (&«/<**)*  (14) 

The  coefficients  c^-i  of  higher  order  are  obtained  direcliy  from  the  equation 

ZotC^lt  - bn  - £ o^cJT’  (15) 

m-t 

whence  the  coefficient  c**  is  adjusted  in  accordance  with  the  transformation 

c<?Uc  <f>  + 2ctc<f2,  (16) 


When  a,  and  6,  are  not  zero,  the  coefficients  cj,1'  of  progressively  increasing  order  are 
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obtained  directly  (rom  the  equation 

£ <w£*»  (17) 

Thue  the  coefficients  are  generated  by  a straightforward  progression. 

The  integration  of  a Taylor  series  in  f is  converted  into  an  integration  of  a power 
series  in  u with  the  aid  of  the  equations 

- -i—  ~ - -- 77  £ (n  + DcSTr1^*  (18) 

du  tn+ldu  m + 1 Hmm 

No  further  multiplication  by  power  series  is  required  for  the  change  of  differential  in 
the  integration.  The  coefficients  of  the  polynomials  which  express  the  powers  of  t in 
terms  of  tht  powei  t of  u are  stored  in  a pair  of  matrices. 


Airy  Function 

The  Airy  function  Ai(z)  is  defined  by  the  equation 


or  by  the  equa'ion 


i r 

Ai(ar)  « — j cos(xf  + Jf*)  dt 


The  first  derivative  Ai(z)  is  given  by  the  equation 


AiW-r,r."' 


and  the  second  derivative  A»"(*)  is  given  by  the  equation 


ct r J 


Addition  of  x to  f*  in  the  integrand  changes  it  into  a function  which  can  be  integrated 
in  finite  terms  as  expressed  by  the  equation 

At  (*)  - — I * 3 dt 

<,H  J.m 

+ ± /-*••»**  (23) 

2n 

The  integration  converges  if  the  limits  of  integration  lie  above  the  real  axis,  in  which 
case  the  Airy  function  is  a solution  of  the  differential  equation 

j <<"(*)  - *At(z)  (24) 
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Let  the  nth  derivative  At**^*)  be  expressed  by  the  equation 

*<•»(*)  «=  £ + £ 6{T‘*Mi'(*)  (25) 

The  coefficients  of  progressively  higher  order  are  generated  by  the  recurrence  equations 

«i’*,-(fc+ (28) 

6<«>  - ai-11  + (k+  DbiV"  (27) 

which  are  started  with  the  single  coefficient  of  zero  order,  a***  -=  1. 

The  expression  of  the  Airy  function  in  terms  of  Bessel  functions  is  derived  on  page 
188  of  Watson's  Theory  of  Bessel  Functions “ 

The  exponential  function  in  the  integrand  of  the  Airy  function  can  be  expanded  in 
an  absolutely  convergent  power  senes  in  *.  If  the  path  of  integration  with  respect  to 
t is  varied  to  pass  inward  toward  the  origin  along  a ray  at  the  phase  angle  In.  and 
outward  from  the  origin  along  a ray  at  the  phase  angle  Jit,  then  the  coefficients  of 
the  power  series  can  be  expressed  in  terms  of  gamma  functions  and  can  be  correlated 
with  the  coefficients  in  the  expansions  of  Bessel  functions  of  imaginary  argument. 
The  Airy  function  is  given  in  terms  of  the  McDonald  function  by  the  equation 

xi 

At(t)  « + ■ Kfi(}**)  (28) 

and  its  derivative  is  given  by  the  equation 

Ai’(s)  • - — ^ Af|(|**)  (29) 

Evaluation  of  the  McDonald  function  is  performed  by  reference  to  a subroutine”’*4  for 
the  complex  Bessel  function. 

The  argument  of  z is  in  the  range  0 to  2n,  whereas  the  argument  for  the  subroutine 
is  in  the  range  -«  to  +tr.  If  the  argument  exceeds  the  range  of  the  subroutine,  then 
the  McDonald  function  is  modified  in  such  a way  as  to  increase  the  argument  by  2n. 
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Figure  1.  Cartesian  Coordinates.  Position  of  a point  P with  respect  to 
the  source  5. 
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K„h  = 1 k0x  = - V 2 «0y  = | *„z  = 0 


figure  4.  Azimuthal  Integration  with  respect  to  Wave  Number  •, 

primary  singularities:  «,  secondary  singularities; , circle  of  radius 

Am,  . boundary  of  convergence;  — . primary  branch,  natural  path; 
— . secondary  branch,  natural  path; . actual  path  of  integration. 


«„(r  + h) 


x = - v 8y 


y = jVx*  + y* 


Figure  6 Computation  Time  alo.ig  the  Critical  Line  of  Wave  Crests.  Curves 

for  k0\  x 2 -1-  y2  equal  to  0(6)32. . trapezoidal  integration;  asymptotic 

approximation;  — . integration  by  parts 
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